function stats=osl_resolve_erf_sign_ambiguity(Sin)

try, tmp=Sin.oat; catch, error('Sin.oat not specified'); end; 
try, tmp=Sin.stats; 
try, first_level_contrasts=Sin.first_level_contrasts; catch, first_level_contrasts=1:size(Sin.stats.cope,3); end; 
 
catch, 
    try,    
        Sin.stats=osl_load_oat_results(Sin.oat,Sin.stats_fname); 
    catch,
        error('Neither Sin.stats nor Sin.stats_fname specified correctly.'); 
    end; 
end;

if(size(Sin.stats.cope,4)>1),
    error('osl_resolve_erf_sign_ambiguity not compatible with TF data');
end;

% build neighbourhoods
blanket=20; %mm
for ii=1:size(Sin.stats.mni_coord,1),
    mc=Sin.stats.mni_coord(ii,:);
        
    % Find the nearest index of the beamformed voxels to the specified mni_coord 
    dists=(sum((Sin.stats.mni_coord-repmat(mc,length(Sin.stats.mni_coord),1)).^2,2)).^(1/2);

    voxinds{ii}=find(dists<blanket);
    %count(ii)=sum(voxinds{ii});
end;

niter=6;
count=zeros(niter,length(first_level_contrasts));
for it=1:niter,
    it
    
    for coni=1:length(first_level_contrasts),        
        
        con=Sin.first_level_contrasts(coni);

        for gc=1:size(Sin.stats.stdcope,5),
            for jj=1:length(voxinds),
                tc_jj=Sin.stats.cope(jj,:,con,1,gc);

                vox_coord=voxinds{jj};
                tcs=Sin.stats.cope(vox_coord,:,con,1,gc);
                [a tc_pca]= pca(tcs,1);

                if mean(tc_pca)<0,
                    tc_pca=-tc_pca;
                end;

                ccs=corrcoef(tc_jj,tc_pca);
                if(ccs(1,2)<0),
                    tc_jj=-tc_jj;
                    count(it,coni)=count(it,coni)+1; 
                    if(0),    
                        if(it>1)
                            figure;plot(normalise(tcs'));ho;plot(normalise(tc_pca),'k','LineWidth',2);plot(normalise(tc_jj),'r','LineWidth',2)
                            title(num2str(Sin.stats.mni_coord(jj,:)))
                            pause;
                            close gcf;
                        end;     
                    end;
                end;

                Sin.stats.cope(jj,:,con,1,gc)=tc_jj;
            end;
        end;
        
        count(it,coni)
    end;
end;

stats=Sin.stats;
stats.count=count;
